# a serious of test 
# 
# Author: guochun
###############################################################################

source("./sourceAll.R")

data=getData("bci5.txt")
availible=data$status=="A" & data$sp == "LONCLA" & !is.na(data$gx) & !is.na(data$gy)
subdata=data[availible,]
ph=getData("PH.dbf")
k=getData("K.dbf")
mg=getData("Mg.dbf")
ph=im(ph,xcol=seq(10,990,by=4),yrow=seq(10,490,by=4))
k=im(k,xcol=seq(10,990,by=4),yrow=seq(10,490,by=4))
mg=im(mg,xcol=seq(10,990,by=4),yrow=seq(10,490,by=4))

covr=list("pH"=ph,"K"=k,"Mg"=mg)
windinf=setwin(data=subdata,covr=covr,plotdim=c(1000,500))
plotdim=windinf[[1]]
subdata=windinf[[2]]
covr=windinf[[3]]

population=new("population",x=subdata$gx,y=subdata$gy,plotdim=plotdim)
AgModelObjects=constructModel(population, covr, select=T, pvalue=T, models=list("Poisson","LGCP","Thomas"))

#for(mo in AgModelObjects){
#	agModeling(mo)
#}

#test LGCP proceduces
mo=AgModelObjects[[2]]
re=agModeling(mo)
del=variableSelection(re,0.05)
if(!is.na(del))
	mo=constructModel(population, re@covr[-del], select=TRUE, pvalue=T, models=list("LGCP"))
re=agModeling(mo)

#debug(backwardSelecter)
re=backwardSelecter(re)

#debug(agregativeResidualTest)
agregativeResidualTest(re,100)
varianceDecomposite(population,covr,select=TRUE, models=list("Thomas"))


#test Thomas proceduces
mo=AgModelObjects[[3]]
re=agModeling(mo)
re=backwardSelecter(re)
varianceDecomposite(population,covr,select=TRUE, models=list("LGCP"))

#test with a simulated LGCP population 
en.filter=covr[[1]]
en.filter$v=covr[[1]]$v*0.05+covr[[2]]$v*0.03-covr[[3]]$v*0.004

sigma2true=2
#expected number of points
N=1000
#get a simulated population generated by LGCP
X=rLGCP(en.filter,sigma2true,100,N,plotdim)

source("./sourceAll.R")
population=new("population",x=X$x,y=X$y,plotdim=plotdim)
re1=varianceDecomposite(population,covr[1:3],select=TRUE, models=list("LGCP"))
tildeZtrue=var(as.numeric(en.filter$v))
vhtrue=tildeZtrue/(tildeZtrue+sigma2true)
vhesti=as.numeric(re1[1]/(re1[1]+re1[3]))



#test with a simulated Thomas population
sigma=50
sim.popu=rThomas(5e-4,sigma,20,win=owin(c(0,plotdim[1]),c(0,plotdim[2])))
en.filter=covr[[1]]
en.filter$v=covr[[1]]$v+covr[[2]]$v*100+covr[[3]]$v*5
en.filter$v=en.filter$v/max(en.filter$v)
sim.popu=rthin(sim.popu,en.filter)
plot(sim.popu)
tildeZ=var(log(as.vector(en.filter$v)))
tildeY=sigma^2

source("./sourceAll.R")
population=new("population",x=sim.popu$x,y=sim.popu$y,plotdim=plotdim)
re1=varianceDecomposite(population,covr,select=TRUE, models=list("LGCP"))

#do DCA first on covariables
library(vegan)
covr.data=data.frame(as.vector(ph$v),as.vector(k$v),as.vector(mg$v))
pca=scores(rda(covr.data),choise=1:3)$sites

covr=covrMaker(pca,size=c(121,246),xcol=seq(10,990,by=4),yrow=seq(10,490,by=4))
windinf=setwin(data=data.frame(gx=population@x,gy=population@y),covr=covr,plotdim=c(1000,500))
plotdim=windinf[[1]]
subdata=windinf[[2]]
covr=windinf[[3]]
population=new("population",x=subdata$gx,y=subdata$gy,plotdim=plotdim)
re2=varianceDecomposite(population,covr,select=TRUE, models=list("LGCP"))


